%% main.m
clear;
close all;
clc;

dt = 1e-3;
nt = 2800;
nlyr = 3;
damp = 0.3;
% nx must be smaller and smaller
nx = [1201, 721, 601];
ny = [601 , 201, 301];
dx = [3, 5, 6];
dy = [3, 5, 6];

rho = [2100, 2100, 2100];
vp  = [2500, 2500, 2500];
vs  = [1000, 1000, 1000];

% nlyr = 1;
% damp = 0.3;
% nx = 1501;
% ny = 601;
% dx = 3;
% dy = 3;
% rho= 2100;
% vp = 2500;
% vs = 1000;


for it = 1:nlyr
    [model(it),layer(it),fields(it)] = ...
        setupParameters(dt,nt,nlyr,nx(it),ny(it),dx(it),dy(it),damp,rho(it),vp(it),vs(it),it);
end

dframe = 50;
nframe = 45;

%% set source
% the layer where single source is set
% source.lyr = 1;
source.lyr = 3;
source.vx  = 0;
source.vy  = 0;
source.sxy = 0;
% source.sxx;
% source.syy;
% src_x = floor(nx/2);
% sigma_x = 1 * dx;
% sigma_y = 1 * dy;
% src_x = 750;
% src_y = 490;
src_x = 400;
src_y = 100;
sigma_x = 6;
sigma_y = 6;

% Ricker STF
f0 = 5;            % 主频(Hz)
t0 = 1.0/f0;       % 延迟时间(s)
t = (0:nt-1)*dt;
src_time = (1 - 2*(pi*f0*(t - t0)).^2) .* exp(-(pi*f0*(t - t0)).^2);

% 生成高斯权重矩阵
[X, Y] = meshgrid(1:model(source.lyr).nx, 1:model(source.lyr).ny);
gauss_weights = exp( -((X' - src_x).^2)/(2*(sigma_x/model(source.lyr).dx)^2) ...
    - ((Y' - src_y).^2)/(2*(sigma_y/model(source.lyr).dy)^2) );
gauss_weights = gauss_weights / max(gauss_weights(:)); % 归一化


% 初始化变量
gif_name = '3block_homo_filter4_lowsrc.gif';
delay_time = 0.1;
is_first_frame = true;

%% 并行化
% if isempty(gcp('nocreate'))
%     parpool('local', nlyr); % 使用与nlyr相等的核心数（需物理核支持）
% end

%%
figure(1);
for it = 1:nt
    current_src = src_time(it) * gauss_weights;
    source.sxx = current_src;
    source.syy = current_src;
    % parfor i = 1:nlyr
    for i = 1:nlyr
        % fields(i) = calc_blk_field(layer(i), model(i), fields(i), source);
        fields(i) = calc_blk_field_chrt(layer(i), model(i), fields(i), source);
        fields(i) = chrt_absorb(fields(i),layer(i));
        % fields(i) = abs_bdry(layer(i), model(i), fields(i));
    end
    for i = 1:nlyr-1
        % disp("chrt bdry");
        [fields(i), fields(i+1)] = ...
            chrt_bdry(fields(i),fields(i+1),model(i),model(i+1),layer(i),layer(i+1));
    end

    downsampled_fields1vx = ...
        imresize(fields(1).vx, [nx(3), floor(ny(1)*nx(3)/nx(1))], 'bilinear');
    downsampled_fields2vx = ...
        imresize(fields(2).vx, [nx(3), floor(ny(2)*nx(3)/nx(2))], 'bilinear');

    split_line = zeros(601,10);
    result_vx = [downsampled_fields1vx, split_line,downsampled_fields2vx, split_line, fields(3).vx];
    % result_vx = fields(1).vx;
    if mod(it,10) == 0
        imagesc(result_vx');
        colorbar;
        xlabel('x');
        ylabel('y');
        clim([-1e-8 1e-8]);
        axis equal;
        axis tight;
        title(['v_x场, 时间步=', num2str(it)]);
        drawnow;

        if mod(it,50) == 0
            frame = getframe(gcf);
            im = frame2im(frame);
            [imind, cm] = rgb2ind(im, 256);
            if (is_first_frame == true)
                imwrite(imind, cm, gif_name, 'gif', ...
                    'Loopcount', inf, 'DelayTime', delay_time);
                is_first_frame = false;
            else
                imwrite(imind, cm, gif_name, 'gif', ...
                    'WriteMode', 'append', 'DelayTime', delay_time);
            end
        end
    end

end